### start necessary libraries
# data import libraries
library(xlsx)
library(foreign)
library(rgdal)
library(tidyverse)
library(dplyr)
# reproducible research library
library(knitr)
# plotting and maps libraries
library(rmapshaper) # for simplifying maps
library(ggplot2)
library(leaflet)
library(ggmap)
### set figure save path to figure/
knitr::opts_chunk$set(fig.path = 'figure/')
knitr::opts_chunk$set(fig.width=12, fig.height=8)
Population data in Australia 2009 to 2014, grouped by Statistical Area Level 3 (SA3) from http://www.health.gov.au/internet/main/publishing.nsf/Content/PHN-Demographic_Data
Primary Health Network (‘PHN’) map data from http://www.health.gov.au/internet/main/publishing.nsf/content/phn-boundaries
Concordance table between PHN geographical area and Statistical Area Level 3 from http://www.health.gov.au/internet/main/publishing.nsf/Content/PHN-Concordances
General Practice (‘GP’, also known as family medicine) workforce by PHN found from Australian Government Health Workforce Data website (http://hwd.health.gov.au/), direct links for General Practitioner data by PHN is not available without login and setting search terms. Pre-downloaded 2016 data.
## OGR data source with driver: ESRI Shapefile
## Source: "PHN_boundaries_AUS_May2017_V7.shp", layer: "PHN_boundaries_AUS_May2017_V7"
## with 31 features
## It has 8 fields
## Integer64 fields read as strings: OBJECTID
Places Population, General Practitioner numbers and Population per General practitioner information into map file
Creates a data frame with
# choose all ages statistic
select.age.population <- select(population.2014, SA3_NAME_2011, AllAges)
# calculate population in each PHN
phn.population <- sa3.phn.concordance %>%
merge(select.age.population, by = 'SA3_NAME_2011') %>%
# merge population numbers by SA3 areas
mutate(PHN_Number = as.numeric(AllAges)*as.numeric(RATIO)) %>%
# calculate population in each SA3 as proportion in each PHN
# some SA3 are spread across more than one PHN
group_by(PHN_NAME) %>%
# PHN can have multiple SA3
summarize(population = sum(PHN_Number))
# add GP information to each PHN
# calculate population per GP ratio in each PHN
combined.oz.phn.description <- oz.phn.description %>%
merge(gp.numbers, by = 'PHN_NAME') %>%
merge(phn.population, by = 'PHN_NAME') %>%
mutate(pop.gp.ratio = population/gp) %>%
# just select the variables that we need
select(PHN_NAME, gp, population, pop.gp.ratio)
oz.phn.map@data <- merge(oz.phn.map@data, combined.oz.phn.description, by = 'PHN_NAME')
‘Click’ on region to reveal Population and GP numbers
oz_loc <- geocode('Australia', source='dsk')
## Information from URL : http://www.datasciencetoolkit.org/maps/api/geocode/json?address=Australia&sensor=false
phn.map <- ms_simplify(oz.phn.map) %>% # reduce points to about 5% of original
leaflet() %>% # the map generator
addTiles() %>% # base map, usually OpenStreetMap
setView(lng = oz_loc$lon, lat = oz_loc$lat, zoom = 4) %>% # set view
addPolygons(color="#444444", weight=1, smoothFactor=0.5,
opacity=0.8, fillOpacity = 0.15,
popup = ~paste(PHN_NAME, '<br>Population :',
sprintf('%.0f',population),
'<br>GPs :', gp,
'<br>Population to GP :',
sprintf('%.0f',pop.gp.ratio)),
fillColor = ~colorQuantile("YlOrRd", pop.gp.ratio)(pop.gp.ratio),
highlightOptions =
highlightOptions(color='green',weight=2,bringToFront = TRUE))
phn.map
‘Click’ on region to reveal Population and GP numbers